Collapse dynamics of trapped Bose-Einstein Condensates 
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We analyze the implosion and subsequent explosion of a trapped condensate after the scattering 
length is switched to a negative value. Our results compare very well qualitatively and fairly well 
quantitatively with the results of recent experiments at JILA. 



o 
o 

(N 
o 

Q 

(N 



> 
O 

(N 



c3 



i 

C 

o 
o 



Bose-Einstein condensates (BECs) in trapped ultra- 
cold atomic gases [|l]-|4j are strongly influenced by the 
atom-atom interactions. These interactions are char- 
acterized by a single parameter, the s-wave scattering 
length a. For a > the interatomic potential is repulsive 
and the condensate is stable. On the other hand, if a < 
the interaction is attractive and a uniform condensate is 
unstable against local collapses. The trapping potential 
stabilizes a condensate with a sufficiently small number 
of particles N < N c , where N c is a critical value below 
which the spacing between the trap levels exceeds the at- 
tractive mean- field interparticle interaction [0 . Trapped 
condensates with a < have been obtained in experi- 
ments with 7 Li at Rice [|[ , and recently at ENS j| . 

Over the last decade, the creation of condensates with 
tunable interparticle interactions (tunable BEC) has at- 
tracted a great deal of interest. There have been several 
proposals on how to modify the scattering length a ■ 
The idea of varying the magnetic field and meeting Fes- 
hbach resonances Q has been successfully implemented 
for Na condensates at MIT ||. However, the Na experi- 
ment was limited by large inelastic losses |l(J. Recently, 
a tunable BEC of 85 Rb atoms has been realized at JILA, 
without large particle losses p|,pT|. These experiments 
constitute an excellent tool for analyzing the influence of 
interatomic interactions on the BEC properties. 

The JILA experiment 11 allows for the creation of 
large condensates with a = ai n u > and a subsequent 
sudden change of the scattering length to a = a co u apse < 
0. After this change of a, the condensate undergoes an 
implosion (collapse) followed by the ejection of relatively 
hot atoms (burst atoms), which form a halo surrounding 
a core of atoms at the trap center (remnant atoms). 

The collapse of a spatially homogeneous condensate is 
described by the well-known nonlinear Schrodinger equa- 
tion (NLSE) which in the context of BECs is called the 
Gross-Pitaevskii (GP) equation. This collapse has a va- 
riety of analogies, such as self-focusing of wave beams in 
nonlinear media, collapse of Langmuir waves, etc. (see 
e.g. |L^-|l4|] and refs. therein). The collapse of the so- 
lutions of NLSE has been extensively investigated and 
it has been found that the dimensionality of the system 
plays a crucial role M JlI . In 3D one has a weak collapse 



where the singularity is reached at a finite time. Before 
this happens the collapsing cloud is described by a uni- 
versal Zakharov solution p5| . This solution consists of 
quasistationary tails and a collapsing central part. The 
number of particles in the collapsing part continuously 
decreases, whereas the density increases. The dynamics 
of collapse in the presence of dissipation has been also 
analyzed (see [ p~3|Jl4| and refs. therein). The dissipation 
is introduced through a nonlinear imaginary (damping) 
term in the NLSE, which prevents the appearance of the 
singularity if the nonlincarity is at least quintic. 

The collapse of a trapped condensate has been recently 
analyzed in several theoretical papers |lr^- |20| , ^| . Kagan 
et al [fnl argued that three-body recombination should 
be explicitly included in the GP equation as an imagi- 
nary loss term. The recombination "burns" only part of 
the condensed atoms and prevents a further collapse of 
the cloud once the peak density becomes so high that 
the three-body loss rate is already comparable with the 
mean-field interaction. Then the collapse turns to ex- 
pansion and the trapped condensate can undergo macro- 
scopic oscillations accompanied by particle losses. Kagan 
et al considered the case of a comparatively large recom- 
bination rate constant, where a single collapse does not 
have an internal structure. By using the formalism of 
Ref. fL7|] , Saito and Ueda (l^] have observed rapid in- 
termittent implosions of the collapsing BEC cloud. This 
resembles the distributed collapse discussed by Vlasov et 
al |]l3f in the context of collapsing cavities in plasmas. 
Saito and Ueda estimated the energy of the atoms re- 
leased during the explosion, and predicted the formation 
of nonlinear patterns in the course of the collapse J2^| . 

A different approach was suggested by Duine and Stoof 
|filf , who proposed binary collisions as the source of burst 
atoms in the JILA experiments. Finally, Kdhler and 
Burnett (22) have recently analyzed the JILA collapsing 
condensates with the help of a non-Markovian nonlinear 
Schrodinger equation, suggesting that the burst atoms 
can be formed due to the violation of the common s-wave 
scattering approximation. 

In the present paper, we analyze the implosion and 
subsequent explosion of the BEC in the conditions of the 
recent experiments with 85 Rb at JILA (ll]]. A detailed 
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analysis of measurable quantities is provided by numeri- 
cal simulations of the GP equation which includes three- 
body recombination losses as proposed in Ref. p7| . Our 
results agree fairly well with the data of JILA. 

We consider a Bose-Einstein condensate of initially N 
atoms of mass m confined in a cylindrically symmetric 
harmonic trap. We restrict ourselves to the trap em- 
ployed in 85 Rb experiments, i.e. a cigar-shaped trap with 
axial frequency uj z — 2n x 6.8Hz, and radial frequency 
Lu p = 2tt x 17.5Hz ||,[llj. Assuming a sufficiently low 
temperature and omitting the presence of an initial ther- 
mal cloud, the behavior of the condensate wavefunction 
-0 is governed by the GP equation (cf. [jl7|l ) : 

m = f-^ + V(r) + g\M 2 ^ T ^-\^A ^ U) 
\ 2m 12 ) 

where V(r) = m{ui 2 p 2 + lo 2 z 2 )/2 is the trapping poten- 
tial, g = 47r?i 2 a/m, and ip is normalized to N. The last 
term in the rhs of Eq.(Q) describes three-body recom- 
bination losses. The quantity L3 is the recombination 
rate constant for an ultra-cold thermal cloud, and the 
numerical factor 1/12 accounts for the reduction of the 
recombination rate by a factor of 6 in the condensate. 

As in the conditions of the JILA experiment, we con- 
sider an initial scattering length ai n u > 0. For this value 
of a we obtain the ground-state condensate wavefunction 
by evolving Eq. (Q) in imaginary time. At t — the scat- 
tering length is abruptly switched to a value a co u apse < 0. 
The simulation of the subsequent dynamics under the 
conditions of JILA involves a very demanding numeri- 
cal procedure, due to very different time and distance 
scales at t = 0, during the collapse, and after the ex- 
plosion. In our simulation we have numerically solved 
Eq.(|l]) by means of the Crank-Nicholson algorithm with 
variable spatial and time steps. We have taken a special 
care of the time and spatial numerical discretizations and 
checked that our results do not change significantly when 
a more accurate sampling is used. 

An additional problem in simulating the BEC dynam- 
ics from Eq.(|l|) is the lack of knowledge of exact val- 
ues of the three-body rate constant L3. The experi- 
ments [p3[ based on the measurement of losses in ther- 
mal clouds set an upper bound for L3, due to difficul- 
ties to distinguish between two- and three-body losses. 
The rate constant L3 is safely determined only far from 
the Feshbach resonance (154. 9G for 85 Rb). A value of 
L 3 = 4.24 x 10~ 24 cm 6 /s has been measured at 250G 
where a — — 336ao p3fl . The existing predictions for L3 
are related to the case of large positive a |2^j2^]. For 
the recombination to deeply bound states, which is the 
case at a < 0, the predictions contain a number of phc- 
nomenological parameters f2q| . In our simulations we 
rely on the experimental value of L3 at 250G and assume 
the dependence L3 oc a 2 to obtain L3 for magnetic fields 
in which \a\ is smaller (a < 0). This particular choice 



is within the error bars of the JILA measurements of L3 
p3| and leads to a fairly good agreement between our 
calculations and the JILA experimental results |ll| for 
the dynamics of collapsing condensates. 

As discussed above, the condensate collapses if N is 
larger than a critical value. In the absence of three-body 
losses, the cloud collapses continuously and the central 
density approaches infinity at a finite time t co u apse . Af- 
ter some time from the start of the collapse, the cloud 
becomes spherical and is described by the universal Za- 
kharov solution |l5| . The size of the central part of the 
cloud decreases as (tcoUapse—t) 1 ^ 2 and the central density 
increases as 1 / (t co u apse —t). In this stage of the collapse 
the presence of the trapping potential is not important 
(see ^0|). We have tested the appearance of the universal 
Zakharov solution in our simulations. 

The presence of three-body recombination changes the 
situation drastically. Once the central density becomes 
such that g\tp(0,t)\ 2 ~ (hL 3 /12)\ip{0, i)| 4 , the collapse 
stops since the recombination losses prevent further in- 
crease of the central density |l7j . We have checked that 
for most of realistic values of L3 the Zakharov solution 
is not realized in the course of the contraction to maxi- 
mum density. For comparatively small L3 the maximum 
density of the cloud is rather high, and the number of 
particles in the central part of the cloud is small. These 
particles are rapidly burned by the recombination and 
the central density drops. However, the central region 
is then quickly refilled by the flux of particles from the 
wings of the spatial distribution. Therefore one obtains 
a set of intermittent collapses, i.e. the collapse becomes 
distributed |jl|,[l!|. As each intermittent collapse burns 
only a very small number of atoms the total number of 
particles presents a smooth time dependence. 

In our calculations we have analyzed the implosion and 
successive explosion of a condensate for the initial num- 
ber of atoms N — 6000 and N — 15000, and for a co u apse 
ranging from — 25ao to — 300ao- We have considered 
o-init — for the case of N = 6000, and a^n — 7ao for 
N = 15000, in order to compare our results with those 
at JILA jyj. Typically, we observe that the conden- 
sate contracts mostly radially and reaches a maximum 
central density after a time t co u apse which ranges from 
0.5 ms to several ms. Then intermittent collapses oc- 
cur. Close to the maximum density in each intermittent 
collapse, the central region of the collapsing condensate 
becomes spherical. As expected, the collapse stops when 
(^3/12)^(0, t)\ 4 ~ g\ip(0, t)\ 2 |l6). At this maximum 
density the total number of condensed atoms (N to t) de- 
creases due to recombination losses. Due to the presence 
of a set of intermittent collapses, the time dependence of 
N tot shows a step- wise decay. However, the average over 
short time intervals of the order of the time interval be- 
tween neighbouring intermittent collapses allows us to fit 
Ntot by an exponential exp(— t/rdecay)- After the BEC 
explodes, we observe the formation of a dilute halo of 
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burst atoms surrounding a central cloud of atoms. This 
reproduces qualitatively the picture observed at JILA. 



FIG. 1. The fractions Nburst/N (circles) and N re mnant/N 
(squares) in % versus a co llapse for N = 6000, a ini t = 0. The 
solid curve shows N cr /N . 
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FIG. 2. The times t C oltapse (upper figure) and Tdecay (bot- 
tom figure) versus a C oiia P se for N = 6000, a init = (circles), 
and N = 15000, ainit = 7&o (squares). In the upper figure the 
triangles show the experimental results of JILA for TV = 6000. 



We calculate various quantities: t co u aps 
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number of burst (N{, U rst) and remnant (N remnant ) atoms, 
and the axial and radial energy of the burst atoms. We 
determine the number of burst atoms by integrating the 
condensate density over the axial coordinate and fit- 
ting the tail of the obtained radial profile by a Gaus- 
sian. The axial and radial energies per particle in the 
burst were calculated by averaging the axial and ra- 
dial Hamiltonians, H z = — 7i 2 V 2 /2m + mw^ 2 /2 and 
Hp = — fi, 2 V 2 /2m + mujpP 2 /2, over the density distri- 
bution. Since the presence of the remnant cloud may 
introduce errors in the determination of the burst ener- 
gies, we have excluded the central region from the average 
of H z and H p . We prevent possible errors by extracting 



central regions of different widths ranging from one to 
several harmonic oscillator lengths. After a typical sim- 
ulation time of 20 ms, we have checked that our results 
for the burst energy per particle are independent of the 
width of the excluded central region. 

In Fig. [l] we present the fractions Nb urst /N and 
N remnant I N versus a coUapse for N = 6000 at a time of 
20 ms. After this time Nburst and N remnan t reach sta- 
tionary values. In the same figure we depict the criti- 
cal value N cr = kaho/\a C oiiapse\, were k — 0.46 is the 
stability coefficient for 85 Rb p7|], and a% = y/h/mc5 
with Q = (WpWz) 1 / 3 . A similar picture has been obtained 
for N = 15000. The burst fraction Nburst/N varies be- 
tween 15% and 25% and only weakly depends on N and 
(^collapse- This is in good agreement with the results of 
JILA PJ, where Nburst/N ~ 20%. One can also see that 
Nremnant > N cr , which is expected and is in agreement 
with the experiments at JILA. 

Fig. | displays the dependence of t co u apse and T decay 
on a co uapse for N = 6000 and N — 15000, respectively. 
As observed, neither characteristic time changes signifi- 
cantly with changing N. The time of collapse ranges from 
0.5 to 4 ms for considered values of a co iiapse- The decay 
time Tdecay weakly depends on a co llapse- For N — 6000 
it ranges from 2.5 ms at a co iiapse = — 25ao to 1.6 ms at 
acollapse = -300a . For the same values of a co Uapse at 
N = 15000, the time Tdecay ranges from 2.2 to 1.1 ms. 
These results are in excellent agreement with Ref. |ll[ . 

Fig. U shows the radial and axial burst energies versus 
^collapse for N — 6000 and N = 15000. Both energies 
smoothly depend on a co u apse and N. The radial energy 
is of the order of 100 nK, and the axial one of the order of 
50 nK. We have performed simulations for a large range of 
values of a and and found that the energy of the burst 
atoms is proportional to a 2 /Ls, as observed by Saito and 
Ueda This is expected as the burst energy should 

be of the order of the maximum mean-field interaction 
g\ip{0, t)\ 2 at the trap center for t ~ t co u apS e, and this 
interaction has exactly the same scaling. 

For the three-body rate constants used in our calcula- 
tions, the results are in a fair agreement with the data 
of the JILA experiment. However, we are not able to 
reproduce the JILA results for the radial burst energy in 
the case of N = 15000 and a co iiapse — 100a . The JILA 
radial burst energy for this particular set of parameters 
is by a factor of 4 larger than that observed in our simu- 
lations, whereas the rest of measured energies depart by 
less than 50% from our results. 

To summarize, we have analyzed the collapse dynam- 
ics of trapped condensates after the scattering length is 
switched to a negative value. Our analysis, based on the 
GP equation with three-body losses explicitly included, 
explains qualitatively and to a large extent quantitatively 
the experiments performed at JILA. 
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Radial energy, N=6000 



axially symmetric trap of JILA (second version of p8[). 
They calculated t co ii apse and found the number of burst 
and remnant atoms as functions of the initial number of 
atoms N at a given value of a co u apse . Their results for 
tcoiiapse agree very well with ours. 
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Radial energy, N= 15000" 
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Axial energy, N= 15000 * 
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FIG. 3. Radial and axial energies versus a co u apst , for 
N = 6000, a in u = 0, and N = 15000, a mit = 7a . Our 
numerical results (circles) are compared with the experimen- 
tal data at JILA (squares with lines indicating the error bars. 

We acknowledge support from Deutsche Forschungs- 
gemeinschaft (SFB 407), the Alexander von Humboldt 
Foundation, the TMR Network 'Coherent Matter Wave 
Interactions', the Dutch Foundations NWO and FOM, 
and from the Russian Foundation for Fundamental Re- 
search. L. S. wishes to thank the ZIP Program of the 
German Government. We acknowledge fruitful discus- 
sions with S. L. Cornish, E. A. Donley, J. L. Roberts, E. 
A. Cornell, C. E. Weiman, and M. Lewenstcin. 

Note added: After this work was completed we learned 
that Saito and Ueda extended their analysis, also based 
on the GP equation with three-body losses, to the case of 
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